function [gbest,gbestval,fitcount] = EGSA_func(fhd,Dimension,Particle_Number,Max_Gen,lowerbound,upperbound,pos,bool_Plot,varargin)

G0 = 2.8;
g = 0.13;
epsilon = 1e-3;
rand('state',sum(100*clock));
me = Max_Gen;
ps = Particle_Number;
D = Dimension;

[pub, plb, vub, vlb] = SetBound(D,ps,lowerbound,upperbound);
% pos = VRmin+(VRmax-VRmin).*rand(ps,D);

fitcount = ps;
vel = zeros(ps,D); % initialize the velocity of the particles
acc = zeros(ps,D); % initialize the acceleration of the particles

fitness = feval(fhd,pos',varargin{:});
sort_fitness = [fitness' (1:ps)'];
sort_fitness = sortrows(sort_fitness, 1);
fitbest = sort_fitness(1,1);
fitworst = sort_fitness(ps,1);
gbestindex = sort_fitness(1,2);
gbest = pos(gbestindex,:);
gbestval = fitbest;

if(bool_Plot)
    figure(3)
    haxes = plot( 0 , 0 );
    XArray = [1];
    YArray = [gbestval];
    title('EGSA');
end

for iter = 2:me
    if(fitbest == fitworst)
        mass = repmat(1/ps,ps,1);
    else
        mass = (fitness'-fitworst)/(fitbest-fitworst);
    end
    mass = mass / sum(mass);

    G = G0 * exp(1 - (iter / me));
    numKbest = ceil(ps * (1 - (iter /me)));
    numKbest = (numKbest>1)*numKbest + (numKbest<=1)*1;
    for n = 1 : ps
        acc(n,:) = zeros(1,D);
        for k = 1 : numKbest
            force = zeros(1,D);
            indexKbest = sort_fitness(k, 2);
            force(1,:) = pos(indexKbest,:) - pos(n,:);
            R = norm(force(1,:));
            force(1,:) = force(1,:) * rand() * mass(indexKbest) / (R + epsilon);
            acc(n,:) = acc(n,:) + force(1,:);
        end
    end
    acc = G * acc;

    mean_pos = zeros(1, D);
    std_pos = zeros(1, D);
    for d = 1 : D
        mean_pos(1,d) = mean(pos(:,d));
        std_pos(1,d) = std(pos(:,d));
    end

    numescape = ceil((ps-numKbest)/10);
    numescape = (numescape>1)*numescape + (numescape<=1)*1;
    for n = 1 : ps
        bool_escape = false;
        if(iter > 1)
            for e = ps : -1 : (ps-numescape+1)
                bool_escape = bool_escape | (n == sort_fitness(e,2));
                if(bool_escape)
                    break;
                end
            end
        end

        if(bool_escape)
%             for d = 1 : D
%                 if(std_pos(1,d) < 1)
%                     pos(n,d) = rand();
%                 else
%                     if(rand() >= 0.5)
%                         pos(n,d) = mean_pos(1,d) + (3 * std_pos(1,d));
%                     else
%                         pos(n,d) = mean_pos(1,d) - (3 * std_pos(1,d));
%                     end
%                 end
%             end
%             if(rand() > 0.5)
                center = zeros(1,D);
                for k = 1 : numKbest
                    indexKbest = sort_fitness(k, 2);
                    center = center + pos(indexKbest,:);
                end
                center = center / numKbest;                
                for d = 1 : D
                    S = rand();
                    S = (S>=0.5)-(S<0.5); % S = (S>=0.5)*1+(S<0.5)*(-1)
                    pos(n,d) = center(1,d) + S * g * abs(center(1,d) - pos(n,d)) * log(1/rand());
                end
%             end
            vel(n,:) = zeros(1,D);
        else        
            vel(n,:) = vel(n,:).*rand(1,D) + acc(n,:);
            vel(n,:) = ((vel(n,:)>=vlb(n,:))&(vel(n,:)<=vub(n,:))).*vel(n,:)+(vel(n,:)<vlb(n,:)).*vlb(n,:)+(vel(n,:)>vub(n,:)).*vub(n,:);
            pos(n,:) = pos(n,:) + vel(n,:);
            pos(n,:) = ((pos(n,:)>=plb(n,:))&(pos(n,:)<=pub(n,:))).*pos(n,:)+(pos(n,:)<plb(n,:)).*plb(n,:)+(pos(n,:)>pub(n,:)).*pub(n,:);
        end
    end

    fitness = feval(fhd,pos',varargin{:});
    sort_fitness = [fitness' (1:ps)'];
    sort_fitness = sortrows(sort_fitness, 1);
    fitbest = sort_fitness(1,1);
    fitworst = sort_fitness(ps,1);
    gbestindex = sort_fitness(1,2);
    gbest = pos(gbestindex,:);
    gbestval = fitbest;
    
    if(bool_Plot)
        XArray = [ XArray iter]; 
        YArray = [ YArray gbestval];
        set( haxes , 'XData' , XArray , 'YData' , YArray );
        drawnow
    end
end
%     if(mod(iter, 10) == 0)
%         iter
%     end
end